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We explore the nematic order of model goethite nanorods in an external magnetic field within 
Onsager-Parsons density functional theory. The goethite rods are represented by monodisperse, 
charged spherocylinders with a permanent magnetic moment along the rod main axis, forcing the 
particles to align parallel to the magnetic field at low field strength. The intrinsic diamagnetic 
susceptibility anisometry of the rods is negative which leads to a preferred perpendicular orientation 
at higher field strength. It is shown that these counteracting effects may give rise to intricate phase 
behavior, including a pronounced stability of biaxial nematic order and the presence of reentrant 
phase transitions and demixing phenomena. The effect of the applied field on the nematic-to-smectic 
transition will also be addressed. 
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I. INTRODUCTION 

The effect of external fields on the phase behavior of 
anisometric colloids (rods and plates) has received con- 
siderable attention over the last years^in particular in the 
domain of density functional theory The presence of 
an interface or wall breaks the translational and orienta- 
tional symmetry and leads to local structures drastically 
different from those found in the bulk |^ |^ ^ |1| . Macro- 
scopically different behavior may be brought about for 
instance by a gravitational field H, 0- If the buoyant 
mass of the colloids is sufficiently high, inhomogeneous 
density profiles are built up along the vertical dimension 
of the system. In some cases, this may lead to the for- 
mation of multiple phase equilibria not encountered at 
zero-field |lll|l3. 

In addition, much effort has been put into investigating 
the behavior of anisometric colloids in an applied elec- 
tric or magnetic field. Owing to the fact that the elec- 
tric polarizability (or magnetic susceptibility) is different 
along the short and long axis of the particle, electric or 
magnetic dipole moments are induced which give rise to 
an additional energetic contribution to the free energy. 
The resulting competition between minimal self-energy 
and maximal configurational entropy of the rods drasti- 
cally changes the orientational structure of the system 
and leads to qualitatively different phase behaviour. An 
important difference with the previously mentioned class 
of fields is that an external electric or magnetic field is 
only coupled with the orientational degrees of freedom of 
the rods or plates and therefore does not directly lead to 
spatial inhomogeneities. 

The first systematic attempt to incorporate the effect 
of these directional external fields into the classic Onsager 
theory (IL] for lyotropic anisometric particles has been 
reported by Khokhlov and Semenov [l^l- Within a sim- 
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plified variational approach, general phase behavior sce- 
narios were presented for both hard and semi-fiexible rod 
systems subjected to various types of directional fields. 
Later on, similar studies have been carried out within a 
numerical treatment of the Onsager theory [l^l fll 'l and 
Lee-Parsons [3| and MWDA-type 16] density functional 
approaches. In the latter case various (spatially) inho- 
mogeneous liquid crystal states which may occur at high 
packing fractions, e.g. smectic and (plastic) solids, were 
also considered. 

Although the experimental study of electro-magnetic 
field effects on colloidal sus pen sions has been pioneered 
a long time ago (see Ref. [13 and references therein), 
the topic has been subject of renewed interest because 
of recent ex peri ments on colloidal goethite (a-FeOOH) 
suspensions |l7l Il8l| . These systems consist of charged, 
bar-shaped nanorods with peculiar magnetic properties. 
These particles not only possess a permanent magnetic 
moment directed along their longitudinal axis, originat- 
ing from uncompensated surface spins within the anti- 
ferromagnetic crystal lattice, but also an enhanced mag- 
netic susceptibility along their short axes. This means 
that additional magnetic moments are induced perpen- 
dicular to the main axis upon applying an external field. 
These unique properties become manifest in particular in 
concentrated, nematic suspensions subjected to magnetic 
fields below 1 Tesla. At low field strengths, the induced 
moments are weak and the permanent ones give rise to 
enhanced nematic alignment along the field direction. 
However, at high field strengths the induced moments 
are dominant and cause the rods to orient with their lon- 
gitudinal axes perpendicular to the field. The associated 
reorientation of the nematic director can be clearly ob- 
served from X-ray scattering measurements |l9j . 

Although the realignment phenomenon can be un- 
derstood directly from the counteracting effects of the 
permanent and induced magnetic dipoles, very little is 
known about the overall phase diagram of goethite sys - 
tems as a function of applied field strength. In Ref. [Tg , 
a first attempt has been made towards a global under- 
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standing of the phase behavior but the analysis there was 
restricted to ideal systems, described by a simple Boltz- 
mann distribution for the rod orientations, and weakly 
correlated systems described by an expansion of the On- 
sager free energy up to first order in the degree of nematic 
order. Both approaches may be used for very dilute 
isotropic suspension where particle interactions are not 
too important, but lack predictive power for dense sys- 
tems where interactions lead to strong deviations from 
the ideal Boltzmann orientation distribution associated 
with the magnetic energy. 

In this paper, we give a full numerical analysis of the 
phase behavior of model goethite suspensions, starting 
from Onsager-Parsons density functional theory. The nu- 
merical approach not only allows us to explore the full 
nematic density range, it also provides a reliable way to 
probe subtle changes in the nematic orientation symme- 
try as a function of the applied field strength, in par- 
ticular transitions between uniaxial and biaxial nematic 
states. We stress that the model goethite system we con- 
sider in this paper is a strongly simplified one. The par- 
ticles are considered as monodispcrse, charged sphero- 
cylinders interacting via electrostatic repulsions. Weak 
attractive van der Waals forces are present in the experi- 
mental systems but are difficult to incorporate the- 
oretically and are not considered here. The interaction 
energy between the total dipole moments on the rods is 
estimated to be of order lO^^fcsT and can therefore 
be safely neglected. Moreover, the particles' consider- 
able size polydispersity (along all three particle axes) not 
only leads to a wide variety of particle shapes, but also 
a strong concomitant spread in the magnetic and elec- 
trostatic properties (e.g. surface charge). Therefore, all 
quantities presented here pertaining to the electrostatic 
and magnetic properties should be considered as typical 
values rather than quantitative averages. 

This paper is constructed as follows. In Sec. II the 
Onsager-Parsons approach will be presented and adapted 
for charged particles and the presence of a directional 
field. In Sec. Ill we quantify the average magnetic and 
electrostatic properties of the goethite rods. Depend- 
ing on the relative contribution of the (average) perma- 
nent and induced dipole moments, several phase diagram 
scenarios for goethite were constructed. They will be 
discussed in Sec. V. In the next section, we scrutinize 
the implications of the magnetic field on the nematic-to- 
smectic transition by means of a first-order bifurcation 
analysis. In the Appendix we supplement our numer- 
ical work with an analytic variational approach based 
upon the Gaussian approximation. This provides us with 
a simple tractable theory for strongly ordered nematic 
states. 



II. ONSAGER-PARSONS THEORY 

The starting point of our analysis is the magnetic en- 
ergy of a single goethite rod which consists of two parts; 



a contribution for the remanent magnetic moment along 
the main rod axis, linear in the magnetic field strength 
B and one representing the induced magnetisation per- 
pendicular to the main axis which depends quadratically 
on B. Following Ref. J^J, the total magnetic energy can 
be written as 

f3Urr,{cos0) = -JBViicosO) + KB'^V2{cos9) (1) 

in terms of the Legendre polynomials Vn with 9 the an- 
gle between the main rod axis and the direction of the 
magnetic field. The quantities J and K, with dimensions 
and (T is Tesla), respectively, are related to the 
(average) remanent dipole moment /i^ and the diamag- 
netic susceptibility anisometry Ax — X\\ ~ Xi- < of the 
nanorods via: 

J = (3fir, K = f3Axvo/3fXo (2) 

with vq the rod volume and the magnetic permeability 
in vacuum. All quantities will be specified in Sec. III-A. 

As a first approximation, the bar-shaped goethite rods 
are modelled as (uniaxial) spherocylinders with equal 
length L and diameter D, bearing a uniform electric sur- 
face charge. Following Onsager ^j, the charged rods 
interact via an effective hard core repulsion, character- 
ized by an effective diameter Dcs > D which depends on 
the charge density on the particle and the ionic strength 
of the solvent. The Onsager-Parsons free energy of the 
system in the presence of an external directional field can 
be cast into the following functional form [Til l2l| : 

+ (/?[/„ (cos 0))^ (3) 

where the brackets denote an orientational average ac- 
cording to some singlet orientation distribution function 
(ODF) /(ri) characterizing the average orientational con- 
figuration of the system in terms of the solid angle Vl. 
Here, c and (p denote the effective (dimensionless) num- 
ber density c = {■k/A)NL'^Dcs/V and packing fraction 
= Nv^s/V with weff = {■K/^)LDlf^ + {■k/&)DIs the ef- 
fective volume of the spherocylinder. The last term in 
Eq. |(2Jl represents the external magnetic contribution 
[cf. Eq. 1^], whereas a and p quantify the orientational 
and packing entropy, respectively, defined by the follow- 
ing angular averages: 

a[f] - (ln47r/(r!))^ 

P[f] - ^{((--7(^^'"')», + ^ + |(^)'}(4) 

where 7 is the angle between two spherocylinders with 
orientations O and fi'. The second term in p[f] arises 
from end-cap contributions to the (electrostatic) repul- 
sion between two short spherocylinders and is strictly 
speaking only valid in the isotropic state 0|. The con- 
tribution rj expresses the so-called 'twisting effect' aris- 
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ing from the orientation-dependent nature of the electro- 
static interaction plll23| : 



7^[f] = -((-Sin7(r!,f2')ln[sin7(f^,f^')]))/ 
-(ln2-l/2)p[/] 



(5) 



The importance of this effect is quantified by a twist- 
ing parameter h — k^^ / D^s, defined as the ratio of the 
Debye screening length and the effective rod diam- 
eter |23|. All properties pertaining to the electrostatic 
interactions will be specified in Sec. III-B for the case 
of goethite. The function gp{4>) originates from the Lee- 
Parsons |25|, rescaling of the original second virial 
theory and is proportional to the Carnahan-Starling ex- 
cess free energy (in units fc^T per particle) for a 
hard sphere fluid |27| via: 



5P(0) 



1 - (3/4)0 
(1-0)2 



(6) 



The equilibrium ODF is determined by applying a formal 
minimization of the free energy. This yields the following 
self-consistency condition: 



f{n) = Zexp 



cgp{<j)) / ujin,n')f{n')dn' 



X exp [— /3C/ni(cos0)] 



(7) 



where Z is obtained from the normalization condition 
/ f{^l)d^l = 1. w is the orientation-dependent part of the 
second virial coefficient for two charged spherocylinders 
at fixed solid angles fl and O' [2ll| : 

x{l-h (ln[sin 7(r2, fl')] + In 2 - 1/2)} (8) 

The solid angle fl is conveniently parametrized in terms 
of a polar angle < 6 < n and an azimuthal one < (/? < 
27r, so that dil = sinOdif. Throughout the remainder 
of this text, 9 always refers to the angle between the 
rod main axis and the direction of the magnetic field. 
The azimuthal angle ip then describes the projection of 
the rod axis onto the plane perpendicular to the field. 
Eq. {TI) is solved iteratively according to a discretization 
scheme outlined in Ref. |23| using a 2D grid of angles 
{di,'Pj}i,j=i,N of mesh size N > 30. 

To specify the orientational symmetry in the various 
nematic states we define the following nematic order pa- 
rameters: 



Sn = (^.(cosfl))^, n=l,2 



A 



^sin^ 9 cos 2iy9) 



/ 



(9) 



The first one, Si, quantifies the degree of dipolar order 
due to the permanent dipoles at finite field strengths. 
Note that for non-dipolar rods 5i = and f{ft) is in- 
variant with respect to the inversion 9 t: — 9. The 
second one, S'2, is usually associated with the nematic 



order parameter and measures the (quadrupolar) orien- 
tational order of a nematic state along the direction of 
the magnetic field. Finally, A is nonzero only if there 
is a preferential direction of alignment within a plane 
perpendicular to the field. In this case, the system pos- 
sesses two mutually perpendicular nematic directors (one 
parallel to the field and one perpendicular to the field) 
and the nematic state is therefore biaxial (BX). If A is 
zero, the rods' projections are distributed randomly in 
the azimuthal plane and the system is of uniaxial (U) 
symmetry. 

The sign of S'2 is also important and allows one to dis- 
tinguish two types of nematic order. First, if < S'2 < 1 
the particles are preferentially aligned along the field, 
corresponding to common 'po/ar' nematic order. Alterna- 
tively, if— 1/2<S2<0 the particles are mainly oriented 
in a plane perpendicular to the field direction, leading 
to 'planar' (or anti- nematic) order. Note that the latter 
type of nematic order only occurs in systems subjected to 
disorientational external fields and is not stable at zero- 
field conditions p^ . 

Once the ODF has been obtained, the thermodynamics 
and phase behavior of the system can be inferred from 
the osmotic pressure and chemical potential. These are 
conveniently expressed in terms of the parameters cr, p 
and 77 [cf. Eqs. H4I5|) ] via: 



ATT , 2 dcpgp 

pli — c + c 



{p[f] + hv[f]} 



(3fl = \nc + a[f] + 2c[gp + ^^] {p[f] + h^[f]} 



+ {f]U^{cos9))j 



(10) 



At phase coexistence, these quantities must be equal in 
each of the coexisting phases. Second order phase tran- 
sitions from e.g. uniaxial to biaxial nematic symmetries 
can be localized by means of a first-order bifurcation ana- 
lysis, as discussed in detail in Ref. p9| . 



III. INTRINSIC ROD PROPERTIES 
A. Magnetic properties 

The electronic and magnetic properties of t he g oethite 
rods have been extensively discussed in Ref. jig- Here, 
we shall only briefly recall some of the basic quantities 
we need as input for our calculations. First of all, the 
remanent magnetic dipole moment of the rods /i,- is es- 
timated to be 10"^ fiB (M-B is the Bohr magneton). The 
diamagnetic susceptibility Ax at room temperature is 
1.7 10"'^ and the average particle volume vq = 5.6 10"^'^ 
m^. Using these numbers in Eq. Q we obtain J = 2.28 
and K = 0.72 T'^. These values need to be con- 
sidered with some care because the magnetic properties 
are size-dependent and the inherent size polydispersity 
of the system therefore leads to a considerable spread in 
J and K. 
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To check whether these numbers are representative for 
the experimental goethite system we may trace 5*2 as a 
function of B for a dilute (i.e. isotropic at zero-field) sys- 
tem and locate its zero-point. Beyond this point, the ne- 
matic order parameter is negative which signifies a grad- 
ual change towards the planar-type nematic order found 
at high field-strengths. Experimentally, the zero-point is 
located at 0.35 T whereas theoretically we find 1.809 T, 
irrespective of the density (at least in the para-nematic 
density range, as will become clear later). The large dis- 
crepancy is attributed to the incertainty in K which de- 
pends sensitively on the particle size. We can achieve 
much closer agreement by doubling this value such that 
J = 2.28 T-^ and K = 1.44 T^^. This gives a zero-point 
at 0.552 T, in reasonable agreement with the experimen- 
tal value. In the actual calculations we shall fix J at 2.28 
T^^ and vary K to verify the sensitivity of the phase 
behavior with respect to a change of the (dia)magnetic 
properties. 

B. Electrostatic properties 
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FIG. 1: Generalized para-nematic-nematic phase diagram for 
rods in external directional fields; "1> > in case of an orien- 
tational field, <1> < for a disorientational field. Paranematic- 
nematic binodals are given by solid lines, the dotted one rep- 
resents a second-order phase transition from the para-uniaxial 
to the biaxial nematic state. 



The double-layer potential around a charged colloid 
(with constant surface charge density) can be determined 
from the Poisson-Boltzmann (PB) equation which, in our 
case, must be solved for a cylindrical geometry. At large 
distances, the electrostatic potential tjj around a cylinder 
with diameter D takes the Debye-Hiickel (DH) form : 



PiPe = TKo{kD/2) 



(11) 



where e is the elementary charge and Kq a modified 
Bessel function. The proportionality constant F depends 
upon the surface charge density of the particles. For 
highly charged particles like goethite, the linearized (DH) 
equation cannot be used to obtain F. Instead, the full 
(non-linear) PB equation must be solved to determine its 
value. Approximate but accurate analytical solutions of 
the PB equation for a cylindrical geometry were obtained 
by Philip and Wooding 'so'l which allow for a straight- 
forward calculation of F by means of a simple iterative 
procedure. Once F has been obtained the effective rod 
diameter can easily be calculated from: 



+ 7£-(l/2) 



D II kQ 

(12) 

where 7s is Eulcr's constant and Q the Bjerrum length. 
Using (Tci w 0.2 C/rri^, _D w 15 nm and ionic strength / sa 
4 -10-2 M we find w 1.5 nm and F w 1.0 -10^ for the 
goethite rods. Eq. (|12|l then gives us Des/D « 1.65. For 
the twisting parameter we thus find h ~ 0.063. These 
results indicate that the effect of twist is expected to 
be rather insignificant. The diameter ratio however is 
quite high so that the effective aspect ratio L/Z^cff of 
the spherocylinder is much smaller than that of the bare 
particle. 



IV. RODS IN DIRECTIONAL EXTERNAL 
FIELDS: GENERAL SCENARIOS 

Before discussing the complex phase behavior of the 
model goethite systems, we shall first present two simple 
general scenarios which occur if rodlike particles are sub- 
jected to an external directional electro-magnetic field. 
The general phase diagram for this case has been pre- 
sented in Fig. 1. This diagram is qualitatively similar to 
the one constructed by Khokhlov and Semenov in Ref. 
[1^ and has been recalculated here for infinitely thin 
hard rods {L/D — > 00, /i =0) based on the free energy 
Eq. For convenience, the external energy Eq. 

has been replaced by 



/3?7ext(cOS0) = -$P2(COS0) 



(13) 



in terms of a general field parameter $ with 9 the an- 
gle between the rod main axis and the field direction. If 
$ > 0, the rods prefer to align along the field, and com- 
mon 'polar' nematic order occurs (indicated by "(-I-)"). 
In this case, the field is referred to as having an 'orienta- 
tional' effect on the rods Note that due to the Boltz- 
mann factor, exp[— /3C/m], in Eq- O, the isotropic state 
no longer exists at finite field-strengths. Instead, dilute 
systems now show weakly aligned para-nematic order, 
indicated by "P" . Both nematic phases are of uniax- 
ial symmetry and the first-order para-nematic-nematic 
coexistence region terminates in a critical point above 
which the system changes from one state to the other in 
a continuous fashion. As to magnetic fields, the present 
scenario may be observed for e.g. rods with a positive 
magnetic susceptibility, Ax > 0, leading to an induced 
moment along the main rod axis. The magnetic field then 
gives rise to liquid crystalline order of the orientational 
quadrupolar type [l^ . Obviously, similar behavior is ex- 
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pected for rods with a permanent magnetic dipole mo- 
ment along their main axis, hke goethite (orientational 
dipolar field). In this case one refers to orientational 
dipolar-type order. 

In the opposite case ($ < 0) the field has a 'disori- 
entational' effect and the rods preferentially orient in a 
plane perpendicular to the field. Both para-nematic and 
nematic phases are now of the planar, or anti-nematic 
type, indicated by "(-)". Moreover, in the concentrated 
nematic phase the rods apparently pack more favorably 
if they attain an additional direction of alignment within 
the plane. The nematic phase thus has a biaxial sym- 
metry. For $ < 0, the first-order para-nematic-nematic 
transition terminates in a tricritical point Beyond 
this point, the transition from one state to the other oc- 
curs by means of a second-order phase transition. The 
present type of ordering may occur if the rods have a 
negative diamagnetic susceptibility anisometry so that 
the induced magnetic moments are perpendicular to their 
main axes (disorientational quadrupolar order). This is 
the case for the induced moments of the goethite rods. 
Similar behavior has been found recently for systems of 
colloidal gibbsite platelets in a magnetic field, where mag- 
netic moments are induced along the short axis of the 
particle [3ll |. 

It is clear from the above that the goethite systems 
are expected to display characteristics from both scenar- 
ios. This will become clear in the next section where we 
shall discuss some explicit phase diagrams for our model 
goethite systems. 



V. PHASE DIAGRAMS FOR GOETHITE 
A. Quadrupolar scenario 

Fig. 2a shows a phase diagram for 'goethite' sphero- 
cylinders. This scenario is similar to the experimental sit- 
uation, judging from the location of the dotted line which 
marks the gradual change from polar to planar-type ne- 
matic order in the dilute regime, discussed in Sec. III-A. 
We refer to this diagram as the 'quadrupolar scenario' 
since the high-field region of the diagram is largely de- 
termined by the induced magnetic moments. Hence the 
appearance is similar to the disorientational quadrupolar 
scenario depicted in Fig. 1. 

At low field strengths, the remanent moments domi- 
nate and the diagram is governed by the orientational ef- 
fect of the field, i.e. the para-nematic and nematic states 
are both uniaxial and the rods are strongly aligned along 
the field direction. However, upon increasing the field 
strength, the degree of polar order will decrease since the 
induced moments (perpendicular to the main rod axis) 
become more pronounced. At some point, the uniax- 
ial nematic state changes to a biaxial one and a first 
order (para-)uniaxial-biaxial nematic coexistence devel- 
ops. The coexistence region eventually narrows down 
towards a tricritical point, beyond which the PU^^^-BX 



transition becomes second-order At very high field 
strengths, the induced moments will completely outweigh 
the remanent ones and force the rods to orient almost 
perfectly in a plane perpendicular to the field. The pla- 
nar PU — BX bifurcation then becomes reminiscent of a 
quasi-2D isotropic-nematic transition [T^ . 

Let us now focus on the field-induced transitions cor- 
responding to the homogeneous (mono-phasic) systems, 
given by the horizontal curves in Fig. 2a. In the dilute 
regime, the ODF changes continuously from polar-type 
(peaked around 6 = ) to planar-type (peaked around 
9 = TT /2 ) and the transition can be roughly localized 
from the condition 5*2 = (dashed curve). Note that the 
corresponding curve does not represent a phase transi- 
tion, it merely localizes a gradual change of signature of 
the ODF. The curve is virtually independent of the pack- 
ing fraction since the ODF in the dilute regime is mainly 
determined by the Boltzmann factor in Eq. 0. 

In the concentrated regime, there is a transition from 
the uniaxial to the biaxial nematic state corresponding 
to a hard bifurcation [s^. This is evidenced by the be- 
havior of the nematic order parameters in Fig. 2b, which 
display a distinct jump just above B — 0.4 T. Note that 
this behavior is quite different from a second-order phase 
transition (or soft bifurcation) where A rises continu- 
ously from zero, rather than jumping to a finite value. 
Physically, the order parameter jump can be associated 
with a sudden reorientation of the main nematic director 
from parallel to the field (at low B) to perpendicular to 
the field (at high B). A similar phenomenon has been 
observed experimental^, albeit at a somewhat lower ap- 
plied field, B ~ 0.2 T [ia. 

It is important to note from Fig. 2b that the rods 
remain sufficiently ordered along their main directors 
throughout the entire field range. By rotating its main 
nematic director perpendicular to the field the system 
is able to sustain the level of polar nematic order with- 
out changing to planar nematic order such as in the di- 
lute regime. This particular property allows us to per- 
form an asymptotic analysis of the free energy, valid for 
strongly ordered polar nematic states. In Appendix A, 
we shall use Gaussian trial ODFs to approximately locate 
the transition and gain some insight into the underlying 
mechanism. In the limiting case of perfectly aligned rods, 
valid for very dense nematic systems, the analysis yields 
a simple expression for the free energy difference between 
the C/^"''' and BX phases at fixed concentration, which 
we will derive in an alternative way here. 

If the spherocylinders are perfectly aligned, the 
phase (denoted by "||") is represented by a collection of 
rods all parallel to the field. For the BX- phase (denoted 
by "_L" ) we assume that all rods point along a nematic di- 
rector perpendicular to the field. At fixed concentration, 
the excluded-volume entropy p (for 7 = 0) is identical 
in both states and the interactions therefore do not con- 
tribute to the free energy difference. Conceptually, the 
system can be considered as an ideal ensemble of dipoles, 
represented by spins. In the ||-state the magnetic energy 
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FIG. 2: (a) Phase diagram of 'goethite' spherocylinders with J = 2.28 r~\ K = 1.44 T"^ and L/DcB = 6.3 {L/D = 10). 
Binodals are indicated by solid lines. The dashed line in the para-nematic region corresponds to S2 = 0. The dotted curve in 
the para-nematic regime denotes a second-order phase transition from the uniaxial P(7'~' to the biaxial BX nematic state, the 
one in the nematic regime represents a hard bifurcation from (/'^' to BX. The dashed curve is the result from the Gaussian 
analysis (Appendix A). The corresponding behavior of the nematic order parameters S2 (solid line) and A (dotted line) is 
shown in (b) for (j> = 0.40. Note the distinct jump at B = 0.442 T. 



of a single spin pUm is then equal to —JB+KB^ (parallel 
to the field) or JB + KB^ (anti-parallel). In the _L-state 
all particles are perpendicular to the field and the mag- 
netic energy is invariant with respect to the direction of 
the spin, hence pUm — —KB^/2 for all particles. The 
free energy difference A!F± ~ -^W i^ ^^'^ easily calculated 
from the spin partition function and reads 

AJ^ = ln[cosh JB] - ^KB^ (14) 

It is easily verified that AT > at low field strength 
whereas AJ- < at high fields, indicating that perpen- 
dicular (_L) order is indeed favored at high field strengths. 
At the realignment transition AJ- is zero, and the corre- 
sponding field strength is found to be i? = 0.513 T for the 
present case. This value represents an upper bound for 
the transition and is not too far away from the numerical 
results. The latter were obtained by numerically mini- 
mizing the free energy difference between the two states. 
The deviations are due to the orientational entropy of 
the rods (neglected in the spin concept) for which we can 
partly account using the Gaussian approximation, dis- 
cussed in Appendix A. 

B. Dipolar scenario 

If we reduce the diamagnetic susceptibility anisome- 
try (by lowering K) the diamagnetic effect becomes rel- 
atively unimportant at low fields. The influence of the 
remanent dipole moments is then expected to govern the 
phase behavior in this regime. We see from Fig. 3a that 
the topology is indeed very similar to the orientational 
scenario of Fig. 1. The para-nematic-nematic coexis- 
tence region terminates in a critical point above which 



the system changes gradually from one state to the other 
without any discontinuity or jump in the associated ne- 
matic order parameters. 

At high field strength, the induced moments become 
dominant and give rise to a reopening of the phase gap 
beyond some critical B-value. Note that the nematic 
phase is now of biaxial symmetry, like in Fig. 2a. A re- 
markable difference with the previous scenario however 
is that the para-nematic phase becomes biaxial as well 
at _B > 2.38 T. In this regime, a coexistence between 
two biaxial phases (a para-biaxial nematic (PBX) and 
a biaxial nematic (BX) one) develops which eventually 
closes off at a critical or consolute point located a.t B = 
2.848 T. Beyond this point the system gradually changes 
from one state to the other, similar to the para-nematic- 
nematic transition above the PJ7^+' — t/^^-' critical point. 
In fact, the para-biaxial-biaxial demixing region can be 
considered as the high-field analog of the PC/'+' — C/^"*"' 
transition. Both involve a coexistence between phases of 
equal symmetry and the entropic mechanism underpin- 
ning the demixing is governed by a competition between 
orientational entropy (favoring the weakly ordered para- 
nematic state) and packing entropy (favoring the nematic 
state). 

An obvious consequence of reducing K is that the tran- 
sitions pertaining to the homogeneous systems shift to 
much higher _B-values, as we see in Fig. 3a. In the con- 
centrated regime, the transition from the uniaxial to the 
biaxial state no longer corresponds to a hard bifurcation 
(or a realignment of the nematic director). Fig. 3b shows 
that the nematic order parameters no longer display a 
real jump but merely a kink at the transition point in- 
dicating that we are dealing with a second-order phase 
transition, similar to the one occurring in the parane- 
matic density regime. 
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FIG. 3: (a) Phase diagram for J = 2.28 T'^ and K = 0.72 T'^ All dotted curves denote second-order phase transitions from 
the uniaxial to the biaxial nematic state, (b) Behavior of the nematic order parameters 5*2 (solid line) and A (dotted line) as 
a function of B for = 0.40. 
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FIG. 4: Phase diagram for J = 2.28 T'^ and K = 1.08 
T-^ All dotted curves denote second-order phase transitions 
from the uniaxial to the biaxial nematic state. A triple line 
is located at B = 1.41 T. 



may be expected upon increasing field strengths. Similar 
to Fig. 3a, the transition from the uniaxial to the biaxial 
nematic state in the nematic density regime corresponds 
to a second-order phase transition, comparable to the one 
depicted in Fig. 3b. 

The present scenario can be nicely connected to the 
previous ones by focussing on the triple equilibrium. 
If we increase K, the concentration of the para-biaxial 
phase (middle dot) is expected to move closer to that 
of the coexisting biaxial phase (right dot) so that the 
biaxial-biaxial region is pushed out of the diagram. At 
some K-vahie both concentrations meet at a critical end- 
point where the PBX-BX region has completely disap- 
peared. From this point on the scenario will be similar 
to Fig. 2a. If we decrease K, the opposite happens: 
the uniaxial-biaxial region is squeezed out at the benefit 
of the biaxial-biaxial region (Fig. 3a). Simultaneously, 
the lower PU'--^^ — binodals detach from the upper 
ones. The latter now constitute a separate 
coexistence region, enclosed by a lower tri-critical point 
and an upper consolute point. 



C. Intermediate scenario 

If we choose an intermediate value for K we get an even 
richer phase diagram as can be seen from Fig. 4. Close 
inspection reveals that the diagram contains features of 
both previous scenarios. This is particularly notable at 
high field strengths, where we observe two para-nematic- 
nematic regions (involving a PU — BX and a PBX ~ BX 
coexistence) reminiscent of the upper regions of Fig. 2a 
and Fig. 3a, respectively. Upon lowering B the two re- 
gions meet at a triple line, indicating a triphasic coexis- 
tence between a uniaxial nematic phase and two biaxial 
nematic phases each with a different concentration. 

Marked reentrant and remixing effects are notable 
around (p = 0.36 where a sequence of phase coexistences 



D. Smectic order: frustration versus realignment 

At packing fractions exceeding roughly 40 % the ne- 
matic phase becomes unstable with respect to smectic or- 
der. To investigate the implications of the magnetic field 
for the nematic-to-smectic transition we have performed 
a stability analysis of the nematic state with respect to 
small spatial density modulations pertaining to smectic 
order. The bifurcation analysis is outlined in Appendix 
B. In Fig. 5 the results are illustrated in terms of a bifur- 
cation diagram for the dipolar scenario, corresponding to 
Fig. 3. At low field strengths, the nematic-smectic bi- 
furcation is virtually the same as without external field, 
as we expect. However, significant changes occur at the 
transition to the biaxial state which, as we see from Fig. 
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FIG. 5: (a) Nematic-smectic bifurcation for the 'dipolar' scenario (J = 2.28 T-\ K = 0.72 T"^). The solid curve corresponds 
to smectic layering along the field (q || B), the dashed branch is associated with layering perpendicular to the field (q 1. B). 
The dotted line marks the second-order uniaxial-biaxial nematic transition. The corresponding normalized layer spacings A,/I/ 
are depicted in (b). The inset shows the evolution of (af) — ((u-ei)^), the root-mean-square projection of the rod vector 
u onto the reference axes {i — x,y,z), with (a^) (dotted), (a^) (dashed) and (a^) (solid). The magnetic field is directed 
along the z-axis. (c) Sketched snapshots representing the nematic phase at the bifurcation density around B = 1.0 T (left) and 
B — 3.0 T (right). The middle one represents the quasi-isotropic structure occurring around B = 1.5 T. 



3b, is associated with a continuous loss of polar nematic 
order (a decrease of 52). As the structure becomes more 
planar-like, the transition to the smectic state persists 
upon increasing field strength but the associated layer 
spacing decreases rapidly and even becomes smaller than 
unity (Fig. 5b). This can be explained from the fact that 
the rods are strongly tilted away from the field direction 
and therefore attain an effective longitudinal dimension 
being smaller than the rod length L. Due to the biaxial 
signature of the nematic reference state, the structure of 
the smectic phase is also expected to be biaxial, although 
a full analysis of the (first-order) nematic to smectic tran- 
sition is probably required to make conclusive statements 
about this. 

At high field strengths, the smectic phase eventually 
realigns in a similar fashion as the nematic phase. In the 
realigned state, the smectic density modulations prop- 
agate along the (main) nematic director perpendicular 
to the field. Similar to the parallel branch, pronounced 
planar-type order occurs (in this case in the xz-plane) 
upon lowering B leading to a decrease of the layer spacing 
and an increase of the bifurcation density. At very high 
field strengths, the nematic-smectic bifurcation becomes 
more or less analogous to the one at low fields. Although 



the reoriented nematic structure remains (slightly) biax- 
ial at high fields, the basic difference between the two 
will be merely a change of the laboratory frame. 

The subtle evolution of the nematic structure upon 
variation of the field strength is perhaps more clearly re- 
flected in the configurational snapshots and the inset of 
Fig 5b, showing the degree of order along the three axes 
of the laboratory frame as a function of B. Note that the 
field points along the z-axis of the frame. At low fields, 
the rods align parallel to field so that (a^) ^ {^x,y)- 
In the reoriented nematic phase at very high fields, the 
rods are mostly pointing along the a;-axis perpendicu- 
lar to the field, hence (a^) ^ (^yz)- There is an in- 
tersection point where (a^) ~ (a^), indicating that all 
rod projections onto the xz-plane have equal probabil- 
ity. The corresponding snapshot shows that the struc- 
ture appears more or less isotropic. However, it is not 
strictly isotropic since the rod vectors only have a small 
y-component, (a^) ^ (a^z)' indicating that they are 
strongly tilted away from the y— axis. Note that in the 
isotropic phase (a^ y z) = l/^i irrespective of the choice 
of laboratory frame. 

Going back to Fig 5, it is not too surprising that the 
two smectic branches meet at exactly the same state 
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point. Both parallel and perpendicular smectic instabil- 
ities are equally favorable because the degree of order is 
the same in both directions. In fact, the quasi-isotropic 
orientation distribution in the a;z-plane can be associ- 
ated with a degeneracy of all smectic instability direc- 
tions within the xz— plane. 

Although the present result may be mathematically 
sound, it remains rather awkward from a physical point 
of view. In particular, one may question whether smectic 
order will be found at all in such weakly aligned config- 
urations. We argue that the present results may be an 
artifact of the assumption that the ODF is constrained at 
the bifurcation point (see Appendix B). The bifurcation 
densities presented in Fig. 4 may therefore deviate from 
those found from the exact eigenvalue equation 33]. In 
the latter case, subtle changes in the ODF are accounted 
for. These may, for instance, cause both branches to ter- 
minate at a terminal point. This particular feature would 
be analogous to the presence of a terminal polydispersity 
for the nematic-to-smectic transition found for parallel 
hard rods with length-polydispersity [s^. Therefore we 
anticipate that the applied magnetic field may give rise 
to a complete frustration of smectic order in a small field 
interval around the intersection point. 

For the intermediate scenario we expect similar fea- 
tures since the nematic structure realigns in an analogous 
fashion compared to the one outlined in Fig. 5. As to 
the quadrupolar scenario, different behavior is expected 
because the sudden reorientation of the nematic direc- 
tor may also take place within the smectic phase. In 
that case, the intermediate quasi-isotropic configurations 
in Fig. 5 are less likely to occur and complete frustra- 
tion of smectic order is probably less pronounced than 
in the other two cases. It would be intriguing to verify 
the possible destruction of smectic order for the various 
cases by means of a computer simulation study. In this 
way more conclusive insight could be generated about the 
implications of an applied magnetic field on the nematic- 
to-smectic transition. 

Although in experiment, smectic order seems to be 
largely suppressed in favor of columnar order (irrespec- 
tive of the apphed field strength H^), smectic textures 
have very recently been observed in between the nematic 
and the columnar phases |3^ . It remains to be investi- 
gated what happens to these textures in an applied mag- 
netic field. 



VI. CONCLUSIONS 

Within the Onsager-Parsons theory we have investi- 
gated the stability of the various nematic states which 
may appear in systems of goethite rods when subjected 
to an external magnetic field. In the present study, the 
goethite rods are represented by charged spherocylinders 
bearing a remanent magnetic moment (leading to pre- 
ferred dipolar order) and a negative diamagnetic sus- 
ceptibility anisometry (leading to preferred planar, or 



quadrupolar, order). These mixed dipolar/quadrupolar 
properties give rise to intricate liquid crystalline phase 
behavior. Depending on the relative contributions of 
the particles' remanent dipole moment and the nega- 
tive magnetic susceptibility anisometry, several scenarios 
were constructed. The quadrupolar scenario, in which 
the effect of the remanent moments is relatively small, 
is characterized by a sudden realignment of the nematic 
director at some critical field strength comparable to ex- 
perimental findings. We have shown that the realignment 
phenomenon can be described appropriately using Gaus- 
sian trial ODFs. 

Upon lowering the susceptibility anisometry, qualita- 
tively different phase diagrams are found. In the dipolar 
scenario, where diamagnetic effects only become mani- 
fest at high field strengths, two separate para-nematic- 
nematic coexistence regions are found at low and high 
fields, the latter involving two biaxial nematic phases. 
At intermediate susceptibilities, a triphasic coexistence is 
found between a uniaxial para-nematic phase and two bi- 
axial nematic ones. Preliminary results for the nematic- 
smectic transition reveal that a similar realignment tran- 
sition may take place for the smectic state. Subtle phe- 
nomena occur in an interval around the realignment field 
strength where smectic order may be suppressed com- 
pletely. 

In the present calculations we have not accounted for 
the bar-shaped geometry of the goethite particles. The 
inherent biaxial shape may have serious consequences for 
the phase diagrams presented here, in particular with 
respect to transitions to the biaxial nematic state. We 
anticipate that biaxial order will be significantly stabi- 
lized because the bar shape makes them prone to biaxial 
nematic order, even at zero- field. Further complications, 
such as size polydispersity could also be addressed by 
considering e.g. binary mixtures of two different-sized 
spherocylinders. Note that the size-dependency of the 
magnetic properties should then also be taken into ac- 
count. However, given the complex phase behavior of the 
monodisperse systems considered here one may question 
whether it is worthwhile to pursue in this direction. 

From an experimental point of view, a promising way 
to reconciliate the present model system with the exper- 
imental one could be to reduce both the intrinsic bar- 
shape and the size-spread of the colloids. The first could 
be achieved by coating the particles with a layer of silica 
which would render them more cylinder-like. The coat- 
ing procedure also opens up the possibility of introducing 
hard-particle interactions by applying a polymer-grafting 
of the silica-coated particles and redispersing them into 
a suitable apolar solvent. However we do not expect this 
modification to give significantly different phase behavior 
since the electrostatic twist effect is of marginal impor- 
tance and all phase diagrams presented here qualitatively 
apply to 'hard' goethite rods as well. The second goal, 
reduction of the polydispersity, can be reached using var- 
ious purification and fractionation methods. In particu- 
lar, reducing the particles' considerable length polydis- 
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persity would be desirable to enhance the stability of 
smectic order. Finally, a systematic variation of parti- 
cle size is expected to influence the relative importance of 
permanent and induced magnetic moments, which would 
offer a means to address different theoretical scenarios. 
These experimental topics are currently under investiga- 
tion at the Van 't Hoff laboratory and significant progress 
has already been made. 



Appendix A : Gaussian approximation 

For strongly aligned states the ODFs are significantly 
peaked around their main nematic directors so that we 
may adopt Gaussian trial functions jl^l to describe the 
orientation distribution in the various nematic states. To 
simplify matters here, we shall neglect the effect of twist 
which, being a minor effect anyway, is expected to be 
even less significant for the strongly aligned systems we 
consider. 

For the uniaxial state (denoted by "||") the following 
Gaussian trial ODF is proposed 



flie) cx exp 



1 



f3Umicos9) 



(15) 



in terms of the variational parameter a. Using Eq. 

in the asymptotic limit (0^1) and some rearranging 

gives: 



fhiO) 



exp [- 



JB 



if < 61 < f 



exp [-^a(-)(7r -0)^ -JB] if f < 61 < tt 

(16) 

where a(±) — a|| ± JB ^ 1 and a\i — a — 3KB^. The 
normalization factor Z is a bit complicated and reads: 




(17) 



With the Gaussian ODF all necessary integrations can 
be performed analytically by means of asymptotic ex- 
pansions for large a\\. After straightforward but lengthy 
calculations we get for the oricntational entropy and mag- 
netic energy, respectively 



with 



(T|[ ~ In a|| — 1 + JB tanh JB — ln[cosh JB] 
+JBa^^'^ (tanh - JBcosh"^ JB) 
= -JBSi + KB^S2 



— Q!ll 

3 II 



(18) 
(19) 



^1 

^2 



tanh JBll- 



- J Ba^^ cosh^'^ J B 



1 - 3a; 



(20) 



containing all contributions up to 0{a^^ ^). For the pack- 
ing entropy p|| we use the leading order contribution of 



the asymptotic expansion performed by Onsager pT| : 

P\\ - 4/\/™^ (21) 

Inserting these expressions into the free energy Eq. Q 
(with h = 0) and minimizing with respect to a\\ yields: 




with 



2JB tanh JB - iKB^ 



(22) 



(23) 



Consistency requires that cgp ^ the regime 

where 5|| < 0. Note that for B — (or cgp ^ 1) the 
regular quadratic dependency q;|| Ac^gp/n is recovered 

A similar treatment can now be given for a biaxial ne- 
matic state (denoted by "_L") which is strongly aligned 
along an in-plane director perpendicular to the field. In- 
troducing as the angle between the remanent dipole 
moment along the particle's main axis and the director, 
the magnetic (Zeeman) energy of the dipole is then pro- 
portional to /I • B = sin ^ sin ip with ip the azimuthal angle 
describing the orientation within the plane perpendicular 
to the director. From Eq. ((T5|l the ODF for the _L-state 
can be written as: 



fci'^'f) exp 



"± ,2 

2 



JBip sin ip 



3KB^ 



(24) 

for ^ <C 1. The biaxial nature of the ODF is evidenced 
by the explicit iy9-dependence. Note that Eq. (|24|l is 
no longer a pure Gaussian (also because of the Zeeman 
term linear in ip) and the angular averages are not easily 
obtained for this case. For large ax we can expand Eq. 
(I24|l to quadratic order in ■(/'■ Applying the normalization 
then yields: 



3^ 

2a± 



1 



3a J 



exp 



2 ^ . 



X 



1 + JB^) sin Lp — 



3ICB^ 



V'^sinV^ (25) 



where K. — K — J^/3. The oricntational entropy in the 
_L-state reads 



a_L ~ Inax - 1 + — ( ^ + ^KB"^ 
a I V 3 2 



and the nematic order parameters are given by 
Si - JBoT^ 



(26) 



^2 ~ (1/2) (3ali - 1) 



(27) 
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up to leading order in a±. The magnetic energy is then 
directly obtained from Eq. H19I) . For the _L-state, how- 
ever, it is more natural to define the nematic order pa- 
rameters along the nematic director fi rather than the 
field direction as in Eq. lf?7l) . This gives: 



5f (1- (3/2)^2)^, 



J a 



1 — 3q! 



(28) 



up to leading order in a±. The latter result shows that 
the biaxial order parameter scales inverse quadratically 
with a± and is therefore very small. Due to the inher- 
ent biaxial structure of the _L-state, the packing entropy 
is difficult to access. Here we assume that the domi- 
nant contribution will be the result for the uniaxial state, 
p± ^ 4:/y^TTa± [cf. Eq. H21|) ] and that biaxiality is ex- 
pected to give higher order corrections of minor impor- 
tance. Minimizing the free energy for the _L-state again 
leads to the form Eq. H22I) but with Si\ replaced by: 



(29) 



Gathering results we arrive at the following free energy 
difference AT — T± — jFy (in units per particle) 

between the two states: 



AT 



ln[cosh JB] - -KB^ + In — 
2 a|| 



4c5p(0) 



-1/2 -1/2 

— q;|I 



— (30) 

a± au 



containing all asymptotic contributions up to 0{a^^). 
For dense systems, cgp{(t>) ^ 1 and a ^ 1 so that the 
first two leading order terms survive [cf. Eq. (|14|l ]. The 
— i3X transition can be localized by solving AT = 
with the aid of Eqs. ^ and (jSHl- The results 

of the Gaussian analysis are shown in Fig. 2a and are 
reasonably close to the numerical results. 



Appendix B : Nematic to smectic-A transition 

To approximately locate phase transitions from the ne- 
matic to the smectic state we may apply a first-order bi- 
furcation analysis starting from the (para)nematic free 
energy Eq. [33Ll3^ . Also here, the effect of twist will 
be ignored, i.e. we set h — Q. A limit of local stability 
for the nematic state with respect to infinitesimally small 
spatial density modulations is then associated with a di- 
vergence of the nematic structure function S'(q). This 
leads to the condition 



the nematic-columnar bifurcation is strongly pre-empted 
by the nematic-smectic bifurcation |39l | and we shall 
therefore not consider the possibility of columnar-type 
instabilities for goethite rods. The shape of the rods 
is enclosed in /m, the cosine-transform of the excluded- 
volume body of two hard spherocylinders, which can be 
deduced analytically [up to 0{D^)] from a straightfor- 
ward but tedious geometric analysis [4Q | . 

We have to realize that in our system smectic layering 
may occur either along the field direction (at small field 
strengths), so that q = {0,0, 1} or perpendicular to the 
field (at large field strengths), in which case q = {1, 0, 0} 
or {0,1,0}. We stress that Eq. l|l?T)) assumes a fixed 
ODF at the bifurcation. At the 'exact' bifurcation point, 
changes in the ODF may also contribute to the loss of 
nematic stability. The corresponding eigenvalue equation 
however is difficult to analyze numerically for biaxial ne- 
matic reference states and we shall not consider it here. 
The bifurcation to the smectic state is given by the wave 
vector q* — 2tt/\* (with A* the characteristic layer spac- 
ing) which gives rise to the smallest physical solution (jf 
ofEq. CT- 

For slender rods {L/D 3> 1) in strongly aligned config- 
urations it is possible to derive simple asymptotic ex- 
pressions for /tv/. Assuming small angular deviations 
from the nematic director, the cosine-transform can be 
written in the the following asymptotic form (henceforth 
D = D,ff)- 

fM{q- n') = -2L^Djl{q/2) sm-f{n, n') - 2nLD^joiq) 

(32) 

with jo{x) = X ^sina; a spherical Bessel function and 
q — qL the rescaled wave vector. Using this result in 
Eq. (|31|l and noting that the angular average of the sine 
is proportional to p[/] we can write [with the aid of Eq. 

my- 



L 



D 



(33) 

where c ~ <j)L/D. At zero-field, a ~ ic^gp/ir and the 
term between brackets reduces to 4. The corresponding 
asymptotic result for the nematic-smectic bifurcation is 
then found at (p* — 0.4037 with corresponding layer spac- 
ing X* / L = 1.293. The bifurcation density is in good 
agreement with the simulation result 0* ~ 0.418 plj . 
We stress that the asymptotic result Eq. H33|l can only 
be applied if the nematic order is strongly polar-like {S2 
close to unity). This is only the case either at very low or 
very high field strengths where the nematic structure is 
not significantly disrupted by an imminent reorientation 
of the nematic director. 



S-\<i) = 1 - 0.gp(0) ((/M(q; n,n')/vo))^ = (31) 

where q is the wave vector pertaining to the direction of 
the density modulations. In principle, these could rep- 
resent smectic or columnar order. For (sphero)cylinders. 
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